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Directed Loop Updates for Quantum Lattice Models 
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Abstract 

This article outlines how the quantum Monte Carlo directed loop update recently introduced can be applied to a wide class 
of quantum lattice models. Several models are considered: spin-s XXZ models with longitudinal and transverse magnetic fields, 
boson models with two-body interactions and ID spinful fermion models. Expressions are given for the parameter regimes were 
very efficient "no-bounce" quantum Monte Carlo algorithms can be found. 

PACS numbers: PACS numbers: 02.70.Ss, 05.10.Ln, 05.30.-d, 75.10.Jm 



I. INTRODUCTION 

The invention of non-local loop updates have made 
quantum Monte Carlo (QMC) simulations an indispens- 
able tool for studying large-scale quantum many-body 
systems. Algorithms with non-local updates are advan- 
tageous to algorithms using only local updates. This is 
because they avoid the low temperature slowing down of 
the configuration selection process that occurs for local 
algorithms which severely limits the accuracy and valid- 
ity of the obtained results. 

The most well-known of the non-local QMC algo- 
rithms is the Loop algorithm |l], ^] which can be used 
directly in continuous imaginary timej3| avoiding the 
Trotter-discretization and has proven efficient for a vari- 
ety of systems. Another method is the Stochastic Series 
Expansion(SSE)Q where one relies on an expansion of 
the exponential in the partition function resembling more 
closely what is done in usual diagrammatic perturbation 
theory. Efficient non-local operator-loop updates in this 
setting was first constructed inQ. A third method is 
the worm algorithm^ first used to measure off-diagonal 
Green functions. This method is very similar to the SSE 
with operator-loops, but the rules described in the origi- 
nal formulation || for moving the worm head differs from 
the rules for constructing the operator-loops. 

Recently it was realized that the rules for constructing 
the SSE operator-loops and the rules for constructing 
updates in the Loop algorithm is in fact just different 
solutions of a set of general equations, the directed loop 
equations, following directly from the requirement of de- 
tailed balance. The particular setting, SSE or space-time 
with continuous imaginary time as used in the Loop al- 
gorithm is in fact irrelevant and a particular solution to 
the directed loop equations can be applied to both cases 
with only minor changes H. Thus the issue is not about 
which method is more efficient -the Loop algorithm or 
the SSE operator-loops. Rather, the issue is how to pick 
the most efficient solution to the directed loop equations. 

In Ref. [pj the directed loop equations of the s = 1/2 
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XXZ model were analyzed in detail. However as men- 
tioned there the directed loop equations apply to a much 
wider class of models and we will here show how to con- 
struct algorithms for general lattice models by giving 
some general solutions to the directed loop equations. 

The directed loop equations possess often many solu- 
tions not all of them giving effective algorithms, and it 
will be important to choose some guidelines for how to 
find effective solutions. In [Q it was emphasized that 
the occurrences of a certain type of move, the "bounce" 
which leads to path back-tracking, effectively undoing an 
update already carried out, should be minimized in or- 
der for the algorithm to be effective. In a region in 
parameter space for the s = 1/2 XXZ model was found 
where bounces can be completely avoided. Here we ex- 
tend this analysis to other models. However also in cases 
where all bounces for a given equation set can be cho- 
sen to vanish, there are for some models still choices to 
be made. In particular this is the case for higher spin 
(s > 1/2) XXZ models[||. Although it is impossible to 
test and compare the efficiency of these choices for gen- 
eral models we have here tested a multitude of choices 
for the s — 1 Heisenberg case. 

Using the solutions of the directed loop equations 
presented here one can construct efficient Monte Carlo 
moves just inputting the matrix elements of the original 
Hamiltonian. This is also the case for the solution of the 
directed loop equations employed in Ref. || , see Ref. || . 
However as we will show, the solutions used here and in 
Ref. lead generally to more effective algorithms. 

To show the versatility of the approach we apply the 
rules here to several systems. Spin-s XXZ models, bosons 
with two-body interactions, spinful ID fcrmions and the 
s = 1/2 XXZ model in a transverse field. 



II. THE LOOP UPDATE 

While it was shown in Ref. [Q that the directed loop 
update applies as well to the Loop algorithm as the SSE 
operator-loop method, we will keep the discussion within 
the SSE formalism here. 

The starting point of the SSE method is the power 
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series expansion of the partition function: 

a n—0 



(i) 



where the trace has been written as a sum over diagonal 
matrix elements in a basis The Hamiltonian is 

written in terms of bond operators H^, where b refers to 
a pair of sites i(b),j(b), 



H = -YH b , 



b=l 



(2) 



where AT 6 is the number of bonds on the lattice. The 
explicit minus sign cancels the minus sign in front of (3 in 
Eq. (Q) and so if all matrix elements of Hi, are positive, 
all terms in Eq. (|l|) are positive. The bond operators are 
further decomposed into two operators; 



Hb = Hi,, 



H 



2,b, 



(3) 



where Hi j, is diagonal and H2,b off-diagonal. 

The powers of H in Eq. (|l]) can be expressed as sums 
of products of the bond operators. Such a product is 
conveniently referred to by an operator-index sequence 



S n = [01,61], [a 2 ,b 2 ], [c 



(4) 



where <Zj S {1,2} corresponds to the type of operator 
(l=diagonal, 2=off-diagonal) and 6, £ {1, . . . , A^} is the 
bond index. Hence, 
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(5) 



where (3 = 1/T. It is useful to define normalized states 
resulting when \a) is propagated by a fraction of the SSE 
operator string: 



\a(p)) ~Y[H aiM \a). 



(6) 



In Fig. ^] a) a particular operator sequence is shown for 
a three-site system. There are three operators, one di- 
agonal and two off-diagonal ones. The state at each site 
is labeled by an integer (encircled) and the propagated 
states can be read of as rows of encircled integers. The 
state \a) is the bottom row. For the discussion of the 
Monte Carlo updates it is convenient to recast the pic- 
ture in Fig. ^ a) into a vertex picture shown in Fig. ^ b) 
where each operator is pictured as a vertex with four legs. 
Each leg is carrying state information and is connected 
to another leg on the same or on another vertex. 

The most important part of a Monte Carlo algorithm 
is how the configuration is updated. In the SSE method 
there are two main types of updates, the diagonal update 
and the loop update. The diagonal update is quite trivial 
and inserts or removes diagonal operators in the operator 
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FIG. 1: a) Operator-sequence for a three site system. There 
is one site for each column in the figure. The states on each 
site are labeled by encircled integers. The operators are shown 
as elongated boxes. There are two off-diagonal operators 
(filled boxes) and one diagonal operator (open box). A prop- 
agated state can be read of as one row of encircled integers. 
\a) is bottom row, while the first propagated state |q(1)} is 
the second row from the bottom, b) the operator sequence 
in a) shown as a vertex picture where all operators have be- 
come vertices, each with four legs. Each leg carry information 
about the state (encircled integer) and is connected through 
a dotted line to a leg on another or the same vertex. The 
dotted lines wrap around the top and bottom of the figure. 
In the vertex picture we do not distinguish between diagonal 
and off-diagonal operators as that is defined uniquely by the 
vertex legs. 



string. It serves the purpose of sampling different lengths 
of the operator string [Q. Here we will be concerned with 
the loop update which changes the type of operators in 
the operator string, but does not change the total number 
of operators. In the Loop algorithm there is no notion of 
the diagonal update and so there the only concern is the 
loop update. 

The algorithm for constructing the loop update is as 
follows. With the configuration mapped onto a linked 
vertex configuration, an initial entrance vertex leg is first 
picked at random. Then the state so on the entrance leg 
is proposed to change into a new state s u with a certain 
probability. An exit- leg on the starting vertex (the vertex 
to which the initial vertex leg belongs) is then chosen 
together with new states for the entrance- and exit legs 
according to a certain probability table. As will be seen 
below it is the solution to the directed loop equations 
which dictates the form of this probability table. The 
probability table is constructed such that the new state 
of the initial entrance leg is required to be equal to the 
proposed state s u . 

Changing the state on the entrance-leg or the exit-leg, 
or both, will result in one or two "link-discontinuities" , 
where states on different legs belonging to the same link 
are different. A configuration with link-discontinuities 
does not contribute to the partition function, so the pro- 
cess must be repeated until the configuration has no more 
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link-discontinuities. 

The process repeats by taking the leg connected to 
the exit-leg of the initial vertex as entrance leg to a new 
vertex, and a new exit leg and state changes are again 
selected according to a probability table. In order not to 
introduce more link-discontinuities the new state of the 
new entrance leg is restricted to be equal to the updated 
state of the previous exit leg. Thus the link-discontinuity 
between the previous exit leg and the current entrance 
leg is removed. A state change of the new exit leg will 
however introduce a new link-discontinuity and so the 
link-discontinuity is effectively moved in front of the path. 

When there is a conservation law such that the state 
change at the exit leg is determined by the state change 
at the entrance leg the link-discontinuities will only van- 
ish when the path closes forming a loop. Then the link- 
discontinuity in front of the path will cancel against the 
discontinuity present on the link on which the initial en- 
trance leg belongs. In contrast, when there is no such 
conservation law a link-discontinuity can vanish just be- 
cause an exit state is not changed, although the entrance 
state was changed. One can then terminate the path 
if there was only one link-discontinuity present before 
this step. This can be achieved by requiring no link- 
discontinuity at the initial entrance leg||. This starting 
condition is not possible when there is a conservation law 
as no new configuration would result, but in the absence 
of a conservation law, state changes on the exit leg can 
occur even if there is no state changes on the entrance 
leg. 

Lets now investigate how detailed balance is satisfied 
for the loop update. We will find the restrictions on 
the probabilities governing the selection of exit-legs and 
states as well as on the initial probabilities. This was 
also done in Ref . |jj , but was restricted to the case where 
there is a conservation law such that the exit state is 
determined by the state change of the entrance leg. In 
general the detailed balance condition reads 



W{s)P{s -> s') = P{s' -> s)W(s') 



(7) 



where W{s) is the weight of the configuration s and 
P(s — > s') is the probability of changing the configu- 
ration from s to s' . For the loop update the probability 
of changing the configuration s = s° to s' = s" can be 
written as a sequence of steps 



P(s^s') = ^i?( S °,e 1 )P s ( S (e 1 ) = So ^s tl ) 

xP(s°,ei 
xP(s\e 2 



s 2 ,x 2 ) 



(8) 
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where R(s°, e{) is the probability for choosing the vertex 
leg ei as the initial entrance leg given the full configu- 
ration s°, P s (s°(ei) = so — ► s u ) is the probability for 
proposing a specific new state s u at the initial entrance 
leg e\. The entrance(exit) leg on vertex i is denoted 



X^ 



X 



x^ 



FIG. 2: A sequence of steps leading to a new spin configu- 
ration. The vertices are shown as horizontal lines while the 
legs are vertical. No states are shown, only the labeling of the 
different legs and the vertices they belong to. 



ej(xj). We denote by s l the full configuration after state 
changes on the i'th vertex in the path, so that s° = s and 
s 1 is the state obtained by possibly changing the states 
at ei and x\. The configuration s n = s' . The notation 
s l (j) refers to the single site state at leg j of the full con- 
figuration s l . P(s t_1 , ei — > s l , Xi) is the probability given 
the configuration s t ~ 1 and the entrance leg ei on vertex 
i to exit the same vertex at xi while changing the en- 
trance state to s l (e,) and the exit state to s l (xi)[[[o|. In 
order not to introduce more link-discontinuities there are 
restrictions on the updated state s l (ei) on the entrance 
leg: s l (ei) = s l ~ 1 (xi^i) for i > 1. For i = 1, the first 
vertex, s 1 (ei) = s u . Thus we might as well substitute 
s 1 (ei) for s u in Eq. (^|). An example illustrating some 
of the symbols used here is shown in Fig. |^. The sum in 
Eq. (||) is over all paths and state changes which lead to 
the new configuration s' . 

The expression for the reverse process where s' is 
changed into s can easily be written down in terms of 
the symbols used in Eq. (|J) as each term in the sum can 
be reversed by just starting at the last exit leg and prop- 
agate backwards until the initial entrance leg is reached. 
Changing the states in opposite order brings back the 
original configuration s. The probability for the reverse 
process can thus be written 



Pis' 



Y,R{s n ,x n )P s {s n {x n ) -» s n -\x n )) 
xP(s n ,x n > s n 1 , e n ) 

XPOS"- 1 ,^-!^"- 2 , 



(9) 



xP(s 1 , xi — > s°, ei). 



Now if one requires detailed balance to hold in the change 
of states on a single vertex 



TF(s m - 1 )P( S m - 1 ,e m 



s m ,x m ) (10) 



= P(s m ,x m ^s m -\e m )W(s m ), 

it follows, by multiplying Eq. (|) with W(s) = W(s°), 
repeated use of Eq. @, and comparison with Eq. (0) 
multiplied by W(s') that detailed balance for the whole 
process is satisfied if 

PWOPsOAei)^ 1 ^)) 
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(11) 

R(s n ,x n )P s (s n (x n ) -> s"- 1 ^)) 

holds in addition to Eq. (|l0|). 

When there is a conservation law, x n and ei refer to 
different legs on the same link and the state changes 
s (ei) — > s (ei) and s n (a; n ) — > s™ (a: n ) occurring in 
Eq. (O) must therefore be opposite to each other. So 
if R(s, ei) is chosen to be uniform independent of s and 
ei, detailed balance requires the probabilities of opposite 
update proposals in the initial step to be the same. With 
no conservation law we should set P s (s°(ei) — ■> s 1 (ei)) — 
8 s or ei Y s i/ ei \, causing no link-discontinuities at the initial 
entrance leg. In this case detailed balance is satisfied 
with a uniform R. 

In addition to the above we should require that the 
path always exits a vertex, 



E E ^- 

%i S'(Xi) 



S , Xi 



(12) 



where the sums are over all possible exit legs Xi and state 
changes on this leg. Note again that s t (ei) is constrained 
to be equal to the exit state s t ~ 1 (xi-\) of the previous 
vertex. When the exit leg equals the entrance leg the 
entrance state is first changed, then the state on the exit 
leg. 

Eq. ( |To| ) involves only the ratio of configuration weights 
for configurations which differ at most by having states 
changed at two legs on a single vertex. (They can also dif- 
fer in the number of link-discontinuities, but these carry 
no weight here.) Because the full configuration weight is 
a product over vertex weights it is sufficient to consider 
each vertex separately. 

To simplify the notation slightly we will hereafter use 
the notation v to mean the state configuration on a single 
vertex. The weight of this single vertex is denoted W(v). 
We also introduce a as 



P(v 



i ,e i -> v l ,Xi) = ^— 



V,Xi) 



wo*- 1 ) 



(13) 



Thus given the values of a it is possible to construct the 
probability tables for how to choose exit legs and exit 
states. The equations governing the values of a, Eqs. (jlC" 
and flU) can be written 

a(v % ~ ,e% — > v\xi) = a(v\xi — > v x ~ ,e») (14) 

]T J2 a ( w<_1 > e « = W -1 ) (1 5 ) 

which constitute the directed loop equations introduced 
in Ref. generalized to the case where there is not nec- 
essarily a conservation law dictating the state change on 
the exit leg. 



III. THE DIRECTED LOOP EQUATIONS 

We will now investigate the structure of the directed 
loop equations, Eqs. (H) and (15). We will first consider 
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FIG. 3: Example of a vertex where the entrance leg is the 
lower left leg. 
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FIG. 4: The vertices (shown below) which results from se- 
lecting the different exit legs (shown above) in Fig. ^. The 
conservation law here is such that the sum of the states below 
each vertex equals the sum above. The + on the entrance 
leg indicates that the state on the inleg is to be increased by 
unity 



situations where there is a conservation law. 

In order to describe the general form of the directed 
loop equations for a general interaction with N\ cgs legs it 
is convenient to abbreviate the labeling somewhat from 
that used in the previous section. To define this new 
labeling, we start by selecting a reference vertex (which 
can be any of the allowed vertices) and label its weight 
W\ . We then choose an entrance leg and label this leg as 
leg 1, and then number the rest of the legs on this vertex 
2, 3 ... n — Megs, see Fig. || for an example of a two-site 
interaction. We then pick a specific way of changing the 
state at the entrance leg. Thus the equations derived 
will apply to the vertex with weight W\ and with the 
specific state change at entrance leg 1. On changing the 
states at both the entrance and exit legs (according to the 
conservation law) one arrives at a new vertex. A specific 
example is shown in Fig. ^. Distributing the weight over 
all possible exit legs according to Eq. dTJ) gives 



Wi = an + au 



ai 7 



(16) 



where we have labeled the weights a^- by their entrance 
(i) and exit (j) legs. Now label the weight of the vertex 
reached by exiting at leg j as Wj. Thus if the exit was 
on leg 2 we would label that vertex Wi. 

Now start with the vertex W<i and change the state on 
leg 2 in the opposite way to what was done when leg 2 
was an exit leg. Exiting on any of the legs, W2 has a 
similar decomposition as W\\ 



W 2 



0-21 + a 2 2 



a 2r . 



(17) 



where now the entrance is on leg 2 on the vertex, with 
weight W 2 , which differs from vertex 1 by having changed 
the states at leg 1 and 2. The weight a%\ corresponds to 
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FIG. 5: The vertices (below) resulting from selecting the dif- 



ferent exit legs on vertex 2 in Fig. 
when the entrance leg is leg 2 and t 
opposite direction to that in Fig. H. 
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le state change is in the 
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FIG. 6: Two different ways of arriving at vertex 3 in Fig. Q 
In the top line the process goes from vertex 1 to 3 via vertex 
2. While in the bottom line it goes directly from 1 to 3. 



the process where the path enters at leg 2 and exits at 
leg 1 . The states are changed in the opposite way to that 
when arriving at W% from Wx, and hence the process is 
undoing the changes and arriving back at Wx- In Fig. ^ 
an example of this is shown. From Eq. ( |l4| ) it follows 
that a2i = ai2. Now one can ask if exiting at leg 3 
or higher yields the same vertex when starting from W% 
as it does starting from W\. The answer to this is yes, 
because starting from W\ one would change the state at 
legs 1 and 3 while starting from W2 one would change the 
states at legs 2 and 3. But W2 differs from Wx only by 
having different states at legs 1 and 2 and thus the state 
at leg 2 is changed twice in opposite directions resulting 
in the same configuration W3. An example illustrating 
this is shown in Fig. ^. The weights are hence uniquely 
defined by this procedure, and one is guaranteed that the 
only vertices which are related by the detailed balance 
equations are those which can be reached by changing 
the state on the entrance leg together with the state on 
any exit leg of the reference vertex. The directed loop 
equations can therefore be written as 



/ an 

ai2 


a.12 ■ 
a-22 ■ 


■ ai n \ 

• a 2n 


1 




/ Wi \ 

w 2 


\ ai n 


0-2n ■ 




w 







(18) 



where the matrix on the left hand side is a real symmetric 
n x n matrix with all entries non-negative to avoid neg- 
ative probabilities. The magnitudes of the diagonal ele- 
ments determine the probabilities for exiting on the same 
leg as the entrance, the so-called bounce processes. The 
bounce processes are generally ineffective as they do not 
change the spin configuration and should be minimized. 



For a given model and parameters there are several such 
sets of equations. 

Although the directed loop equations consist of n — 
Ni e g S equations this number is for many models in prac- 
tice often reduced as one or more of the vertices arrived 
at might not be allowed, that is they have zero weight. 
If so, all entries in the matrix involving transitions to the 
disallowed vertex must be zero, which means that the cor- 
responding column in the matrix is zero. The symmetry 
of the matrix implies that also the corresponding row is 
zero and so the dimensionality of the matrix is effectively 
reduced. An example of this is the s — 1/2 XXZ model 
where only three exit possibilities for each entrance leg is 
allowed, thus in this case the directed loop equation sets 
have dimensionality 3. 

There is a general solution to the directed loop 
equations due to Sandvik|| employed in a number of 
works pT|. This solution which we will label solution A 
reads 



WiWi 



Wx+W 2 



(19) 



which implies that the probabilities take the heat-bath 
form 



Wi W 1 + W 2 + --- + W n 



(20) 



The advantage of this solution is that it is general and 
easy to apply. However it doesn't treat bounce processes 
different from other updates and is thus expected not to 
be as effective in many cases. In Ref. § it was shown that 
there are other solutions which performs much better for 
the s = 1/2 XXZ model. 

The system of equations ( |l8| ) contains n(n + l)/2 un- 
knowns and n equations. When n > 2 there are always 
more unknowns than equations; making many solutions 
possible. However if one seek solutions where all bounces 
are absent, all diagonal elements zero, the number of un- 
knowns is reduced to n(n — l)/2 while the number of 
equations remains n. So without bounces one finds that 
for n — 3 the solution is unique while for n > 3 there are 
again many solutions. While this counting of unknowns 
and equations gives information about when one can ex- 
pect solutions, it does not ensure that the solutions are 
positive which is required in order to have positive prob- 
abilities. 

To investigate the solutions more closely let us start 
with the simplest case where the matrix is reduced to a 
2 x 2-matrix. 



an ai2 

«12 0,22 



W x 
W 2 



(21) 



It is clear that a bounce-free solution an = 022 = can 
only occur when W\ — W2 for which a\2 — W\ implying 
that P(l -> 2) = P(2 -> 1) = 1. When Wx ^ W 2 
bounces are necessary. There are several possibilities for 
writing down a solution in this case, but one solution 
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which also generalizes to bigger matrices is to choose to 
bounce off only the vertex with the biggest weight. If we 
assume W% > W 2 this solution is an = W\ — W2, a 22 = 
and a u = W 2 . This gives P(l -> 2) = W 2 /W u P(2 -> 
1) = 1 and P(l -> 1) = 1 - W 2 /Wi. Note that this 
bounce solution reduces to the bounce- free solution in the 
case W\ = W 2 , thus whenever the directed loop equations 
reduces to a 2 x 2-system one can use the bounce solution 
above. 

For a 3 x 3 system the bounce-free solution (an = 
»22 = 033 = 0) is unique and reads 



012 = {W 1 +W 2 -W 3 )/2 

013 = (W 1 -W 2 + W 3 )/2 
a 23 = (-W 1 + W 2 + W 3 )/2 



(22) 



It is clear that this solution ceases to be a valid solution 
whenever one of the weights is bigger than the sum of 
the two smaller weights. In this case one needs again to 
include bounces. Again this can be done as in the 2x2 
case by only bouncing off the vertex with the biggest 
weight. Assuming that W\ is the biggest weight this 
bounce solution can be summarized as 



an = Wi — W 2 
O12 = W 2l 
ai 3 = W 3 , 



W 3 



(23) 



with all other a's zero. This solution is complimentary 
to the bounce-free solution Eq. (|2|) being valid in the 
regime where the bounce-free solution is not valid. Fur- 
thermore these solutions are continuous at the boundary 
where Wi = W 2 + W 3 . It is also interesting to see how 
these solutions reduces to the solutions described above 
for the 2 x 2-case when the smallest weight (assumed here 
to be W 3 ) goes to zero. Then clearly the bounce solution 
reduces to the 2 x 2-case. At first sight the bounce-free 
solution does not, however its regime of validity shrinks 
when W 3 — ► as Wi > W 2 and the bounce-free solution 
is only valid when Wi < W 2 + W 3 . Thus the bounce- free 
solution is again only valid when W\ = W 2 and so also 
the bounce-free solution reduces to the solution found in 
the 2 x 2-case. 

Going on to the 4 x 4-case it is clear that again we can 
write down the bounce solution by bouncing off only the 
vertex with the biggest weight. In fact this solution can 
be generalized to any n x n-matrix and reads when we 
assume that Wi is the biggest weight 



/ an 

W 2 



W 2 
















/ 



(24) 



where an =W\- (W 2 + W 3 -\ h W n ). This solution is 

valid when one weight is bigger than the sum of the rest 
of the weights. This means in practice that bounces are 
only needed in parameter regimes where one term in the 



Hamiltonian dominates. The probability tables following 
from this solution has a quite simple interpretation. The 
probability for moving between vertices other than the 
one with the biggest weight is zero while that of moving 
from the largest weight configuration to the smaller ones 
is the ratio of the smaller weight to the larger weight and 
unity for the reverse process. The bounce probability is 
unity minus the probabilities for moving to the smaller 
weight configurations. 

Whenever the bounce solution above ceases to be valid 
one can write down a bounce-free solution. However, 
this is not unique as there are four equations with six 
unknowns. The different solutions can be parametrized 
as follows 



«12 
ai3 
a 23 

di4 



(Wi + W 2 - 
(W 1 -w 2 + 
(-Wi + w 2 
W 4 - (a 34 + 



W 3 - W 4 )/2 + a 34 
W 3 -W 4 )/2 + a 24 
■f W 3 + W 4 )/2 - (a 34 
024) 



(25) 



024) 



Here we have assumed that W\ > W 2 > W 3 > W 4 . 
When W 4 — > the requirement that 014 be non-negative 
forces a 3 4 and 024 both to zero. This solution with a 34 = 
«24 = goes smoothly into the solutions for the 3x3 
set described above in the limit W4 — > 0. For later use 
we will term this solution where 034 = 024 = (for all 
W4) for solution B. Solution B has the advantage that 
it is valid whenever the bounce solution above is not. It 
also goes continuously into the bounce solution at W\ = 
W 2 + W 3 + W 4 . 

Solution B implies that the region where one can write 
down a solution without bounces is given by 



Wi + W 2 + W 3 + W 4 > 0, 



(26) 



where W\ > W 2 > W 3 > W4. While this was inferred 
from a special solution this result is in fact general. One 
can show on general grounds that the criterion allowing 
for a bounce-free solution is (for an n x n-matrix) 



W x + W 2 + . . . + W n > 0. 



(27) 



There are many bounce-free solutions when n > 4. A 
general one which reduces to the solution B above when 



W 5 = W 6 
and is 



W n 



is termed solution Bl here 
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(Wi + W 2 ~W 3 ~ W 4 )/2 
(Wi - W 2 + W 3 - W 4 )/2 
(-Wi + W 2 + W 3 + W 4 )/2 
W 4 - W5/2 
(W s - W 6 )/2 

(W„_! - W n )/2 

W n /2 

W 5 /2 



(28) 



G 



a 56 = W 6 / 2 



IV. SPIN-S XXZ-MODELS 



a n -i, n = W n /2 

The requirement that Bl reduces to solution B above 
restricts -W x + W 2 + W 3 + W± > 0. 

The above matrix framework also holds in the case 
where there is no conservation law which dictates the 
state change on the exit leg once the state change on 
the entrance leg is given. Then the matrix dimension 
is increased to n = N\ cgs x Abates, where N\ egs is the 
number of legs on the vertex and Abates are the number 
of allowed state changes on each exit leg. This is because 
one must here also take into account all possible state 
changes on the exit legs. This includes the possibility 
where the state on the exit leg remains unchanged, which 
as described on the previous section will terminate the 
path. The numbers 1 to n thus each have a vertex type, 
a leg and an update type associated with them. Here 
also the state change of the exit leg j in the process ay 
is opposite to the state change on the entrance leg j in 
the process djk- 



Before we study the efficiency of different solutions let 
us consider an example. Consider the spin-s XXZ model 
with nearest neighbor interactions. We have also added a 
magnetic field, an interaction proportional to (S z ) 2 and 
a physically unimportant constant C 

<ij> 

+ J2{-hS* + vSf} 

i 

We take the exchange coupling to be ferromagnetic and 
choose units such that its magnitude is unity. Using the 
normalization of the ladder operators S± = S x ± iS y 



S±\s,m) = y/(sTm)(s±m + l)\s,m±l) (30) 

where m G [— s, s] is the S z quantum number which labels 
the state on a vertex-leg, it is easy to see that the different 
vertex weights for this model arc 



W (n ± 1, to =F 1, n, to) = — y/ (s =p n) (s ± n + l)(s ± m)(s =p to + 1) 
W (n, to, n, to) = C — (j z nm — h (n + m) + v (n 2 + m 2 )^j , 



where the arguments of W represent the states on legs 
1,2,3 and 4 labeled as in Fig. ||, and a = a/Z, where Z is 
the coordination number of the lattice. C must be cho- 
sen such that all diagonal weights are positive. Generally 
we will choose C to be slightly above its minimal value. 
The antiferromagnetic case can be studied on bipartite 
lattices where the spins on one of the sublattices are ro- 
tated an angle 7r about the z-axis giving a sign change to 
the spin-flip-term but not to the J z term. 

For general s the maximal dimensionality of the di- 
rected loop equation sets for this model is 4. However, 
they do not all have dimension 4. For the vertices con- 
taining a spin with the maximum (to = s) spin the equa- 
tion set for updates which attempts to increase this spin 
further has dimension 3 as this attempt will lead to a 
vertex with zero weight. This holds also for updates at- 
tempting to decrease a spin below its minimum value. It 
thus follows that all sets of directed loop equations for 
the s = 1/2 XXZ model has at most dimension 30. For 
s = 1 the update of the vertex shown in Fig. is described 
by an equation set with dimensionality 4. This update 
together with the update where the spin on the same ver- 
tex is decreased, their symmetry related updates (enter- 
ing from another in-leg) and their reverse updates(going 
backwards) are in fact the only updates for the s = 1 
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FIG. 7: Update of a vertex for the spin-1 XXZ model which 
results in the 4 allowed vertices shown below. Here we have 
labeled the single-site states with the eigenstates of S z : m g 
[—1, 0, 1]. From left to right the resulting vertex weights are 
C, C + h — v, 1 and 1. If instead one had considered the 
update where the entrance spin is decreased by one unit, the 
spins on the resulting vertices would have the opposite sign 
and the only change in the weights would be that the second 
vertex from the left would have weight C — h — v. 



model governed by an equation set with dimension 4. 
The other equation sets have dimensions 2 and 3. 

It is interesting to ask for the region in parameter space 
where the ineffective bounce processes can be avoided. 
The defining feature of this region is that Eq. (|2^) should 
hold for all equation sets in the model. For s = 1/2 this 



7 



region was shown in Ref. jTj to be the the region de- 
fined by | J z \ +2|ft,| < 1. That is for XY-like anisotropics 
and moderate fields. To generalize this region to ar- 
bitrary spin-s one can consider increasing the state on 
the lower left leg by unity on a diagonal vertex with 
weight W(n,m,n,rri). The directed loop equations then 
relates this vertex to the vertices W(n + l,m,n+ l,m), 
W(n + l,m — l,n, to) and W(n + 1, to, n, m + 1). In the 
case to = the off-diagonal weights are equal, while for 
to = ±s one of the off-diagonal weights vanishes. Now 
choose C such that the weights of the diagonal vertices 
are always larger than the off-diagonal ones. While this is 
a matter of choice in the SSE method it is always true in 
the Loop algorithm as the diagonal vertices are of order 
unity while the off-diagonal ones are of order the Trotter 
spacing. Then the inequality ( p6| ) takes the form 

| J z m - h + v(2n + 1)| < ^(s + n + l)(s - n) 

x (y/ (s — to + l)(s + m) 

+ v/(s + to+1)(s-to)) 

(31) 

where n G [—s,s — 1] and m G [— s,s]. The most re- 
strictive case is for n — s — 1 and to = s for which the 
bounce-free criterion is 

\J z \s + \h\ + \v\(2s-l)<s. (32) 

The same inequality is obtained by considering lowering 
a spin on a diagonal vertex and when a spin on an off- 
diagonal vertex is changed in the same direction as it 
is being changed by the off-diagonal operator. However 
when a spin on an off-diagonal vertex is changed in the 
direction opposite to how it is changed by the operator a 
stricter condition is obtained, see Fig. ||. This is because 
there are no operators changing the spin on a site by two 
units in the XXZ-model. To avoid bounces in this case 

(s + to) (s — to + 1) = (s + to + 1) (s — to) (33) 

where — s + 1 < to < s — 1. Thus this condition does not 
apply for s = 1/2. The condition Eq. d33| ) is equivalent 
to to = — m, which only is satisfied for to = 0. Thus 
for s < 1 the condition ([33]) does not constrain the no- 
bounce parameter region. However for s > 1 one always 
need bounces for these kind of vertices, and so for s > 1 it 
is not possible to find an algorithm within the framework 
presented here which is entirely without bounces 

For s > 1/2 it is also possible to consider updates that 
change the quantum numbers by more than one unit. 
However because there are no non-zero off-diagonal terms 
where the magnetization at a site is changed by two units 
such an update will be described by an equation set of di- 
mensionality 2. Thus for almost all parameters will this 
update contain bounces and because the path without 
bounces is deterministic (it is determined by the straight- 
through process for diagonal vertices and the diagonal 



n-1 m+1 n m+1 n-1 m+1 n-1 m+2 
n m n+1 m n+1 m-1 n+1 m 

FIG. 8: The vertices resulting from increasing the spin on the 
lower left leg of the off-diagonal vertex, shown left. The right 
two vertices are not allowed in the XXZ model, and so the 
directed loop equation set reduces to having dimensionality 
2. The circles around the leg states are omitted for clarity. 



process for off-diagonal ones) there will be a sizable prob- 
ability for the process to bounce back and forth along the 
predetermined path before it possibly ends by retracing 
its path all the way to the starting point without having 
done any changes. It is thus not expected that inclusion 
of these updates will make the simulation more efficient. 

To demonstrate that different spin magnitudes can 
be simulated efficiently using the same basic code with 
changing just the number of different states on a site and 
the matrix elements of the Hamiltonian we show in Fig. ^ 
the magnetization curves for different 100 sites Heisen- 
berg antiferromagnetic spin-chains with s = 1/2, 1, 3/2 
and 2 at inverse temperature (3J = 100. At low fields 
one can clearly see that the integer spin chains have a 
gap while the half-integer ones are gap-less. At low fields 
the stair-case finite size effects are also clearly seen for 
the half-integer chains. A typical point on the s = 1/2 
curve is based on an average of 10 bins each with 10 4 
MCS (5 x 10 4 MCS for equilibration) and took about 20 
minutes on a single processor 868 MHz Intel Pentium III. 
For comparison a typical point on the s = 1, s = 3/2 and 
s = 2 curves using the same number of equilibration and 
measurement steps took about 2, 5 and 8 hours respec- 
tively. The s — 1/2, 1, 3/2 and 2 simulations involves 6, 
17, 34 and 57 non-zero vertices respectively. 



V. EFFICIENCY 

In Ref. J?j a number of examples for the s = 1/2 XXZ 
model showed that algorithms minimizing the number of 
bounces is generally more effective than those where no 
such minimization is attempted. 

In the s = 1/2 XXZ-model the directed loop equa- 
tions have dimensions 3. Thus in the regime where a 
bounce-free solution exists, it is unique. However, when 
the directed loop equations have dimension 4 or greater 
the bounce-free solution is not unique. The solutions to 
the bounce-free 4-dimensional case can be parametrized 
as shown in Eqs. (^5|). We will now investigate how to 
choose the most efficient of these. This cannot be done in 
complete generality as a general physical model involves 
many equation sets relating the different vertices, and as 
we will see the overall efficiency depends on how the sets 
are interconnected. However we can pick a specific model 
and parameters, and try out different solutions there. For 
this we pick a simple model with a 4-dimensional equa- 
tion set, the s — 1 Heisenberg model, and measure auto- 
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FIG. 9: Magnetization per site and spin magnitude as func- 
tion of h/s for 100 sites Heisenberg antiferromagnetic chains 
having different spins (f3 = 100). The inset shows the behav- 
ior at low fields. The dotted lines are guides to the eye. 



correlations for the staggered magnetization as a function 
of different bounce-free solutions. 

We use the integrated autocorrelation time as a mea- 
sure of efficiency [y. It is a measure for how many sub- 
sequent Monte Carlo steps (MCS) are needed in order 
to obtain statistically independent configurations. A low 
value indicates an efficient algorithm. The integrated au- 
tocorrelation time for a quantity O is defined as 



oo 



t=l 



where Ao(t) is 



AoQb) = 



(34) 



(35) 



where i and t are Monte Carlo times, measured in units 
of one MCS. A MCS is defined here by adjusting the 
number of loop updates during equilibration such that 
every vertex in the linked list is on average visited twice 
(not counting bounces) in one MCS. 

There are only two inequivalent 4-dimensional equa- 
tion sets for the s = 1 Heisenberg model. See Fig. [7 for 
the vertices belonging to one of these sets. The vertices 
belonging to the other set are labeled A', B', C and D', 
and are obtained by changing the sign on all states in the 
corresponding unprimed vertices in Fig.^. 



In the zero field Heisenberg case (J 2 = 1, h = v = 0) 
all weights of the vertices shown in Fig. can be chosen 
to be equal (C = 1). When two or more weights are 
equal there is an ambiguity in the ordering of vertices 
when ordered according to their size as is assumed in 
the specification of the solution in Eqs. (25). However, 
a specific ordering of the vertices and values of 024 and 
034 can be shown to be equivalent to a solution for a 
different ordering of the equal-weights vertices with other 
values of 024 and 034. The explicit relations are shown in 
Appendix A]. Therefore it suffices to choose one ordering 
which we here choose to be ABCD, meaning W\ = Wa, 
W 2 = Wb, W 3 = Wc and W4 = Wd where the letters 
refer to the letters in Fig. 0, and investigate efficiency, 
as measured by the integrated autocorrelation function 
for the staggered magnetization, as functions of 024 and 
034. Because of the time- reversal symmetry we choose 
the ordering A'B'C'D' for the other set and use equal 
values of CI24 and 034 for the two sets. 

To avoid searching the whole two-dimensional space of 
these two quantities we will restrict ourselves to three 
lines where 024 = 1234 = a/2, 024 = a, 034 = and 
«24 = 0, 034 = a respectively. The quantity a is restricted 
to the set [0, W A ] where W A = (-W 1 +W 2 + W 3 + W 4 )/2. 
This restriction follows from the requirement of having 
non-negative weights. 

In Fig. [To] we show integrated autocorrelation times 
for the staggered magnetization for values of 024 and 034 
along the three lines specified above. In the top panel 
all weights in both 4-dimcnsional equation sets are equal 
and it is clearly seen that having 024 = gives the lowest 
autocorrelation times, being close to the minimal value 
0.5 in the whole range of 0,34. Naively one would think 
that for equal weights the most effective solution should 
be the most symmetric one: a.y = 1/3, for all i j. The 
arrow in the top panel shows the autocorrelation time for 
this case. This is clearly not the most efficient solution. 
One should keep in mind that the effectiveness as mea- 
sured here depends in general not just on the rules for the 
directed-loop equation sets seen individually, but also on 
how different sets are are interconnected. The asymme- 
try observed above should probably be attributed to the 
diagonal vertex with an up and a down spin which has 
double the weight of the other vertices. Although it is 
not a part of the two 4-dimensional equation sets it still 
plays a role making an asymmetry in how to choose the 
most effective solution. 

Letting C > 1 all four weights are no longer equal. 
However the two diagonal ones are still equal to each 
other as well as are the two off-diagonal ones. Fig. [l(] 
bottom panel shows autocorrelation times when C = 1.5 
still with the same ordering of vertices as in the top panel. 
Again we see that 024 = is favorable for the efficiency. 
In Ref. Q the spin-1 Heisenberg model with disorder 
was studied using directed loops using the solution where 
all diagonal updates were excluded and other processes 
have equal probabilities. In the language used here their 
solution corresponds to the case where 024 = and 0,34 = 
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1/2 with the orderings ABCD and A'B'C'D'. 

In a magnetic field the time-reversal symmetry is bro- 
ken and the two 4-dimcnsional equation sets are no longer 
equal. Furthermore the diagonal vertices in each of these 
sets are not equal and it suffices to choose the ordering of 
the two off-diagonal ones. Fig. [l| shows integrated auto- 
correlations for the staggered magnetization at h = 0.1 
for different values of 024 and 034. In the upper panel the 
ordering of the vertices are BACD and A'B'C'D' while it 
is BACD and A'B'D'C' in the lower panel. In both pan- 
els we have used the same values of 024 and 034 for the two 
sets. The most efficient algorithm using the ordering in 
the upper panel is found for 024 = and 034 — Wa- For 
the ordering used in the lower panel, the efficiency is seen 
to depend solely on the sum 024 + £134 and is maximized 
for 024 + 034 = Wa- A similar analysis was carried out 
by Harada and KawashimaQ. In the language used here 
their parametrization corresponds to a parametrization 
of the ordering used in the lower panel where 034 = 0, 
and ci24 was varied. The algorithm performing best at 
a field h = 0.1, their solution 1, corresponds to the case 
a24 = Wa which indeed also is among the most efficient 
algorithms found here. 

It is interesting to compare these autocorrelation times 
with those obtained using solution A. The integrated au- 
tocorrelation times for the staggered magnetization us- 
ing solution A corresponding to the top, bottom panels 
of Fig.[l(] and to Fig.ll] are 11.6, 12.5 and 12.0 respec- 
tively, which is significantly more than the most efficient 
algorithms which have autocorrelation times close to 0.5. 

For higher magnetic fields the variation of autocorre- 
lation times with the parameter a is similar to what is 
shown for h = 0.1 in Fig. [ll]. However at high fields the 
differences between the least and most efficient choices 
are less pronounced than in small fields. 
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FIG. 10: Integrated autocorrelation times for the staggered 
magnetization of a 64-site Heisenberg (J z = 1) spin-1 chain 
at (3 — 16 as functions of different values of 024 and 034. The 
order of the vertices when having equal weights is ABCD. 
Solid lines indicate the parametrization a24 = «34 = a/2, 
dashed lines: d24 = 0, (134 = a and dotted lines are 0:24 = a, 
034 = 0. The top panel is for the case where all weights 
are equal, C = 1, h = 0. The bottom panel is for C = 1.5, 
h — 0. The arrow indicates the location of the fully symmetric 
algorithm where a,ij — 1/3 for i 7^ j. 



u 



H — (n(n - 1) + n'(n' - 1)) 



VI. INTERACTING BOSONS 

The directed loops can also be applied to the boson 
Hamiltonian 

H = ^ |— t (a\cij + h.c.^j + vriirij — c\ 
<ij> 



u 



+ -^ni(rn - 1) 



(36) 



which describes interacting bosons hopping on a lattice, 
aj, di and rii are the boson creation, annihilation and 
density operator on site i respectively. The single-site 
states are labeled by the boson occupation number. The 
vertex weights are 



W(n, n', n + 1, n' - 1) = ty/{n + l)n' (1 - S n<n J 

W{n,n',n- l,n' + l) = t^n(n' + 1) (1 - <V,0 (37) 
W(n,n ,n, n) = C — [vnn' — ft, (n + n) 



where a = a/Z, and n,n' G [0, n x ]. C must be chosen 
such that the diagonal vertices are non-negative. How- 
ever, this can not always be achieved as in principle the 
boson occupation number can be arbitrarily large. So in 
practice when U > —\v\ one must put an upper limit, n x , 
to how many bosons can occupy a site. The (1 — S n ,n x ) 
factors are needed to implement this. They prevent the 
boson occupation from increasing above n x . Imposing 
this restriction the boson system has as many states on a 
single site as a spin-s model where 2s + 1 = n x + 1. Fur- 
thermore there is a one-to-one correspondence between 
the non-zero vertices in the boson model and the non- 
zero vertices in the spin-s XXZ model, although their val- 
ues are of course different. Thus the only change needed 
to the spin-s program code is to change the numerical 
values for the weights. A special example of this is the 
well-known mapping in the hard-core limit (U/t — > 00) 
where one can set n x = 1 without penalty and thus the 
boson model is mapped onto a spin-1/2 model. The error 
introduced by imposing a restriction on the boson occu- 
pation number can be controlled by repeating runs with 
different n x 's. One expects only small variations in the 
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FIG. 11: Same as Fig. |l(J except that h — 0.1. In the upper 
panel the ordering of vertices is BACD and A'B'C'D' while 
it is BACD and A'B'D'C' in the lower panel. The arrow 
indicates the best algorithm found in Ref. M. 



results as long as n x is bigger than the typical occupation 
number on each site. In the simulations described below 
we use n x — 4 for which there are 57 non-zero vertices. 

To show the efficiency of solution B applied to this bo- 
son system we show in Fig. [l^ the reentrant behavior of 
the Mott insulating phase in ID first predicted from a 
DMRG calculation (jit). In Fig. |l|we show how the den- 
sity p, compressibility k = dp/ dp = N[3((p 2 ) — (p) 2 ) and 
superfluid density p s varies for a 2D Bose- Hubbard model 
as the chemical potential is increased keeping t/U fixed. 
One can clearly see the insulating regions as well as the 
singularities in the compressibility at the first and second 
superfluid-insulator transition. Because of the relatively 
big value of t/U, being near the tip of the lowest Mott 
lobes, the insulating regions become rapidly narrower, 
and the gaps are relatively small making it necessary 
to lower the temperature further to see the singularities 
at the higher-/! superfluid-insulator transitions. Fig. [l4| 
compares integrated autocorrelation times for the den- 
sity as measured using solution A and B for a smaller 
but similar system to that shown in Fig. [l3|. The differ- 
ences in autocorrelation times are most pronounced at 
low values of the chemical potential and away from the 
peaks seen in the Mott insulating regions which arises 
from the small denominator in Eq. (p5|) . 

To find the regime where it is possible to use bounce- 
free algorithms for this Hamiltonian we consider loops 
changing the number of bosons by unity. The conserva- 
tion law here is that the Hamiltonian conserves the total 
number of bosons, n\ + n 2 = n-j + 77,4. Now increase 




FIG. 12: Density p, compressibility n — dp/dp and super- 
fluid density p s for a N = 128 site ID boson Hubbard model 
with U — 1, p/U = 0.17 as functions of the hopping t/U. 
The inverse temperature is (3U = 320. The inset shows the 
compressibility in the region of the reentrant insulating phase. 
The lower curve in the inset is for N = 256 and f3U = 480. 
The superfluid density shown is measured as the square of the 
boson world-lines spatial winding number^HJ. 

the number of bosons on the lower left leg of a diagonal 
vertex by unity, then the vertices related by the directed 
loop equations are shown in Fig.|l5| and their weights are 
from left to right 

Wi = W(n,m,n,m) 

W 2 = W(n + l,m,n + l,m) (38) 

W3 = W(n + l,m—l,n,m) 

W4 = W(n+l,m,n,m+l) 

where W(ni, 712, 713, 714) refers to Eq. (p7|), and n £ 
[0, n x — 1], m E [0,71a-]. The difference between the diag- 
onal weights is W\ — W2 — vm + Un — p. Now choose 
C such that the biggest diagonal weight (W\ or W2) is 
always bigger than any of the off- diag onal weights W3 
and W4. One then finds from Eq. (Eq) that in order to 
avoid bounces 

\vm + Un — p\ < ty/n + 1 (\fm + Vm + 1 (1 — S m n )) 

where n £ [0,n x — 1], m £ [0, tiJ. The same inequal- 
ity is also obtained when considering lowering the boson 
occupation value on a diagonal vertex. In the hard-core 
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FIG. 13: Density p, compressibility k = dp/ dp, and superfluid 
density p 3 for the 32 x 32 boson Hubbard model with t/U — 
0.03 and fit = 2.5 as functions of the chemical potential p/U. 




FIG. 14: Integrated autocorrelation times for the density p in 
a 12 x 12-Bose-Hubbard model. U/t = 33.333 and fit = 0.5 
for the two solutions A and B. 



n m n+1 m n m n m+1 
n m n+1 m n+1 m-1 n+1 m 

FIG. 15: Increasing the occupation on the lower left leg results 
in the following four vertices. (The circles on the legs have been 
omitted.) 



case where n x — 1 this criterion reduces to |/2| < t and 
\ v ~ A I < t- Changing the boson occupation on a leg on 
an off-diagonal vertex in the opposite direction to how 
the state is changed by the operator results in an addi- 
tional criterion as was the case for spin-s models. Here 
this is 



ty/n = tyn + 1 



(40) 



for n x — 1 > n > 1. Thus this does not apply in the 
hard-core boson case. For n x > 1 the equality is never 
satisfied, and so bounces are always needed for these kind 
of vertices when n x > 1. 



VII. ID FERMIONS WITH SPIN 

The rules described here applies also to ID spinful 
fermions jl6| . However for (anti)periodic boundary con- 
ditions the simulation is restricted to a system where the 
number of spin-ups and downs are both an (even)odd 
number. Thus with these boundary conditions one must 
carry out the simulation at so low temperatures that only 
configurations having a fixed (even)odd number of up and 
down spins contributes. For open boundary conditions 
there are no such restrictions. 

Consider the ID fcrmion Hamiltonian 



H 



E ( c i*Cj* + h - c ) + V(m - l)(nj - 1) 

<ij> 

-Jx (S*S? + SfS]) + J Z S*S* - C) (41) 



where c- a , c. 



and the creation, annihilation and 

density operators of a fermion with spin a on site i. rii = 
n-fj + Tin and Si = c\aci/2 are the particle density and 
spin operator on site i. The states on each site are labeled 
by the charge and spin \q,s) so that |0) = |0, 0),| |) = 
|l,i), I |) = and| U) = |2, 0). The weights of the 

diagonal vertices can be read directly off the Hamiltonian 
above and are 



J 



W(qisx,q2S2,qis 1 ,q 2 s 2 ) = C - i ■ J z sxs 2 + V(qi - 1)(q2 - 1) 



(42) 
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"2 [{qi 



i) 2 + (92 - l) 2 - l] - fl, (si + s 2 ) - A (?i + ga) 



where again a = a/Z and Z is the coordination number. 
As the states will be represented in terms of occupation 
number states, which are bosonic states, one needs a con- 
vention for how the phases are related to the fermionic 
states. We follow ref. (ll| and define 

cJ T |....n iT =0...) = (-l) n n|....n,., T = l...) (43) 
cj i |....n ii =0...) = (-ir + <i\....n td = l...) (44) 

where nf^ is the operator counting the number of parti- 
cles with spin up on sites less than i, and nj is the total 
number operator for spin-up particles. A similar defini- 
tion holds for the annihilation operators. The explicit nj 
for the spin down states is to ensure anticommutation be- 
tween fermi operators with different spin indeces. With 
(anti)periodic boundary conditions 

4+i.t = (-!) Pc M ( 45 ) 

= (~ 1 ) Pc u ( 46 ) 

where P = (—1)0 it follows that for the off-diagonal 
weights to be non-negative nj,, raj must both be 
(even)odd. 

The weight of all off-diagonal vertices where one parti- 
cle, and only one, is transferred from one site to the other 
is t. The off-diagonal vertex where two-particles are inter- 
changed, the spin-flip vertex, has weight VF(jj) = J±/2. 
In all there are 34 allowed vertices (32 if J± =0). We 
consider four types of updates: Adding/removing a spin- 
up/down particle. The conservation law here is the con- 
servation of charge and spin. That is: q\ + q% = + q± 
and si + s 2 = S3 + S4. 

This conservation law is in fact so strong that no 4- 
vertex relations exists. This can be seen by considering 
the process of adding a spin-up particle to the lower left 
leg on a diagonal vertex. In this case the update where 
the exit is on the upper left leg, the "continue straight 
through", will be allowed as well as the bounce. In ad- 
dition there is one and only one off-diagonal vertex that 
arises: If the right legs contains an up-spin, then the 
lower right leg is the allowed exit. The upper right leg 
is the allowed exit if the right legs does not contain an 
up spin. Thus only three vertices are related by the di- 
rected loop equations. It is quite clear that this also holds 
for all the other update processes on a diagonal vertex. 
One must also consider adding a spin-up particle to the 
off-diagonal spin-flip vertex. This results again in three 
allowed vertices, but now they are all off-diagonal, two 
with weight t and one with weight Jj_/2. Finally adding 
a spin-up particle to a leg on an off-diagonal vertex with 
weight t results in one of the two above mentioned equa- 
tion sets or in a 2-dimensional equation set with equal 



weights. Thus for this model all the directed loop equa- 
tions have dimensions 3 or less. 

To find the regime where bounces can be avoided we 
employ Eq. ( |26| ) with W4 = 0. Choosing C big enough so 
that one of the diagonal vertices always have the biggest 
weight it follows that 

t> \J z \/4+\V\ + \U\/2+\H z \/2+\fl\ (47) 

in order to avoid bounce-solutions to the equation sets 
where two diagonal and one off-diagonal vertex are re- 
lated. C cancels as it occurs on both sides of the in- 
equality. For the set where three off-diagonal vertices 
are related one must have 

t > Jl/4 (48) 

to avoid bounces ii t < Jj_/2, otherwise it suffices that 
J± is non-negative. 

To compare the efficiency of solution B with solution 
A we have measured the integrated autocorrelation times 
for the charge density wave (CDW) order parameter at 
q = 7r defined as 

Ocuw(q) = (Pc(q)Pc(-q)) (49) 

where 

Pc (q) = (l/A)^e^(n Tr + n ir ) (50) 

r 

as functions of V for U/t = 2, N = 16 and f3t = 32 for 
both solution A and B. While the autocorrelation times 
are quite long in this case considering the relatively small 
system size the solution B is more efficient than solution 
A. The difference in autocorrelation times are largest for 
small values of V. 

Lets now consider another update type where two par- 
ticles are added. Even though the update types consid- 
ered above are sufficient for ergodicity, this update type 
is necessary if one wishes to measure superconducting 
correlation functions (local pairs). Then starting from 
a diagonal vertex it is clear that one cannot reach an 
off-diagonal vertex. Thus there are only two possibili- 
ties, continuing straight through or bouncing. To avoid 
bouncing the weights of the resulting vertices must be 
equal, which means V = [i = 0. Otherwise bounces are 
necessary. Entering an off-diagonal vertex also leads to 
only two exit possibilities, but there both have weight t 
and so the bounce probability can be set to zero. 

VIII. SPIN-1/2 FERROMAGNET IN A TRANS- 
VERSE FIELD 

As an example of a model without a conservation law 
we consider the spin-1/2 XX Z ferromagnet in a trans- 
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FIG. 16: Integrated autocorrelation times for the CDW order- 
parameter for the two solution A and B in a N = 16 sites sys- 
tem at fit — 32, U/t = 2, fj, = J± = J z = H z = as functions 
of V. The inset shows Ocdw(t) measured for system sizes 
N = 16, 32, 64 at low temperatures fit = IN. 



verse magnetic field, that is a field along the x-direction. 

« = - E {( 5 " 5 i + ^ J ^ + c ) 

+h x J2Sf (51) 



The exchange coupling in the spin XY-plane is restricted 
to be ferromagnetic as a 7r-rotation on one sublattice 
which was necessary to obtain an antiferromagnet in the 
zero field case cannot be employed here without introduc- 
ing a minus sign coming from the /i^-term. Alternatively 
one can view this as an antiferromagnet in a staggered 
magnetic field. 

The transverse field introduces vertices where the sum 
of the spins on the lower two legs is not equal to the 
sum on the upper two legs, see Fig. |l7|. This reflects the 
fact that Sf is not a good quantum number in the 
presence of a transverse field. Thus the conservation law 
utilized in Sec. IV cannot be used and we must include 
the possibilities of a state change on just one leg, the 
entrance or the exit leg, keeping the state on the other 
legs unchanged. The path construction in the absence of 
a conservation law always start by keeping the entrance 
leg on the first vertex unchanged, while it stops when the 
state on the exit leg is not changed. 



© 



FIG. 17: Examples of vertices which do not conserve the total 
spin in the Z-direction. In the transverse field s — 1/2 XXZ 
model all such vertices have weight h x /2. (h x = h x /Z). 
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FIG. 18: The possible exit-legs and exit state changes when 
the entrance leg is the lower left leg on the vertex with all 
spins up and the entrance leg is unchanged. The dashed line 
indicates that the state should remain unchanged, while the 
solid line indicates a spin flip. 



To find the directed loop equations here we first look 
at the vertex where all spins are up and the entrance leg 
is unchanged. Then we have the possibilities of exiting at 
one of the four legs as well as changing or not changing 
the state of the exit leg. In all there are eight possibil- 
ities, see Fig. [l8|, thus the directed loop equations have 
dimension 8. While the dimension is rather big there are 
only two different vertex weig hts, W{\\) and W{\\) (the 
vertices with one leg different from the others are all de- 
generate). Next one should take the same vertex, but 
now change the state on the entrance leg. Now there are 
just seven possibilities as changing both the lower legs 
leads to a zero-weight vertex. This procedure should be 
repeated for all vertices, entrance-legs and update types. 

It is interesting to note that the region where one can 
avoid bounces is in fact bigger in this representation us- 
ing a transverse field than found in Sec. IV. Take the 
situation where the lower left leg on a diagonal vertex 
is changed. The vertices belonging to the equation set 
generated in this way are the vertices that would be gen- 
erated without the transverse field plus the four vertices 
with just the in-leg changed and no change of state on the 
out-leg. The no-bounce-criterion is then (assuming that 
C is chosen such that the diagonal vertices are biggest) 



1 h x 
< - +4— . 

2~2 2 



(52) 



One must also consider the sets of the type shown in 
Fig. |is| where the in-leg is not changed. As there in these 
cases are just two different weights out of a total of eight 
the inequality ( |27| ) is always satisfied. 

As explained in Sec. IV there is an ambiguity in solu- 
tion B when many vertex weights are equal. Here we have 
used the convention that when two weights are equal the 
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one with the lowest exit leg comes first. In the rare case 
where the exit legs are also equal, the one with the no 
state change on the exit leg comes first. We have tried 
solution A, Bl and a solution B2 which is particular for 
systems of dimensionality less than or equal to 8: 

012 = (Wi + W 2 - W 3 - Wa)/2 

013 = (Wi - W 2 + W 3 - W 4 )/2 
a 23 = (-W 1 + W 2 +W 3 + W 4 )/2 

014 = W A -(W 5 + W 6 + W 7 + W s )/4 

C115 = «45 = W 5 /i 

am = a 4s = W n /4: 

a 56 = (W 5 + W 6 - W 7 - W 8 )/4 

a 5 7 = (W 5 - W 6 + W 7 - W 8 )/A 

a 58 = W s /2 

a 6 7 = {-W b + W 6 + W 7 + W 8 )/A 



(53) 



(54) 



This solution reduces also to solution B when W5 = Wq = 
. . . = W„ — and is restricted to the region where — W5+ 
W 6 + W 7 + W 8 > 0. 

We found that solutions A and B2 perform better than 
solution Bl in all cases studied. This can perhaps be 
explained by the fact that in Bl many processes vanish 
when the smaller weights are equal which is the case for 
small to intermediate fields. The difference in efficiency 
between solution A and B2 is small in the cases studied 
here. We have been unable so far to find a bounce-free 
solution which clearly outperforms solution A. 

For h x — the 2d XY-model (J z — 0) exhibits a phase- 
transition of the Kosterlitz-Thouless type|l7| as a func- 
tion of temperature where the helicity modulus as mea- 
sured by the second derivative of the free energy to a 
twist 9 in the boundary conditions shows an emerging 
discontinuity with increasing system size at T c . In zero 
field this quantity is efficiently measured as the fluctua- 
tions in the spatial winding number of the loops. With 
a magnetic field term one can still measure the second 
derivative of the free energy with respect to a twist in 
the boundary condition as fluctuations in a "winding" 
number provided one redefines the winding number to 
include a term coming from the magnetic field term in 
the Hamiltonian: 



d 2 F 



T 
~2 



((wl) + (w 2 hy ) 



(55) 



where we have symmetrized the expression in x 
The modified "winding" number Wh„, where a = x,y 



and y. 
is 



N, 



+8 TT 



(56) 
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FIG. 19: The second derivative of the change in free energy 
resulting from a twist in the boundary conditions are shown 
(top panel) as functions of inverse temperature for different 
square lattice system sizes L x L ranging from L = 8 to L = 
1024. The bottom panel shows the plot of d 2 F/d6 2 L which 
collapses onto one curve for all L as well as the magnetization 
obtained using the biggest system size (L=1024). 



where the sum is over all vertices p in the linked list. The 
Kronecker (^-functions contribute whenever the vertex p 
equals the indicated vertex. a(p) is the a-coordinate of 
the site whose state is changed in the vertex p. The 
symbol 5p^ means that the bond to which the vertex 
is attached must lie in the cr-direction, and N a is the 
number of sites in the cr-direction. 

In Fig. |l^ we have measured d 2 F/d8 2 at h x = 1 as a 
function of (inverse) temperature for L x L square lat- 
tices where L ranges from 8 to 1024. For strong fields the 
spins are predominantly in the x-direction, and it is ex- 
pected that a twist in the boundary condition will cause 
d 2 F/d9 2 to depend linearly on L. In the bottom panel 
we confirm this by showing how the scaled curves(scaling 
factor 1/L) collapse onto one curve for all L's. In this 
plot we also show the magnetization for the biggest sys- 
tem size. 



IX. DISCUSSION 

We have shown how to construct probability tables for 
the directed loop update in the SSE method in a model- 
independent way. The construction involves solutions of 
the directed loop equations. These equations arise from 
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considering vertices that can be reached from a certain 
vertex when states on two legs of that vertex are changed. 
For a given model there are many such equation sets. To 
cope with the many sets it is best to generate them in the 
set-up part of the Monte Carlo program. Efficiency is not 
an issue in this construction process as the same proba- 
bility tables will be employed throughout the simulation. 
To construct the probability tables the program should 
go through every entrance leg on every vertex type occur- 
ring in the model under consideration. For each of these 
entrance legs the program also have to consider every 
possible update, as there should be one probability table 
associated with every entrance leg, vertex and update 
type. For each of these entrance legs and update types 
one finds all possible exit legs and exit states, and thus 
the related vertices. The number of vertices with non- 
zero weights reached in this way determines the effective 
dimensionality of the directed loop equations. Having 
found this dimension and weights for the different ver- 
tices one can then pick the probabilities of moving from 



vertex i to j as P(i — ► j) 



j/Wi, where is gotten 



from the general solutions described in Sec. III. 

The solutions are naturally divided into two classes; 
those with and those without bounces: processes where 
the loop back-tracks along its path. Bounces should gen- 
erally be avoided as they are inefficient. The precise crite- 
rion for when bounces can be avoided is given in Eq. (27). 
When this criterion is not fulfilled, bounces are necessary. 
A general solution which then always can be applied is 
given in Eq. (]24|) where there is only one bounce, namely 
bouncing off the vertex with the biggest weight. 

Whenever bounces can be avoided it is likely that there 
exists a bounce-free algorithm which is more effective. 
This is supported by the results injTj and in the exam- 
ples considered here with the exception of the transverse 
field XXZ model where we have not been able to find a 
bounce-free solution which clearly is more efficient than 
solution A. However the author believes that such a so- 
lution exists. 

For directed loop equation sets with dimension < 3, as 
is the case for all the equation sets in the s = 1/2 XXZ 
model as well as for the ID spinful fermion model consid- 
ered here, the bounce-free solution is unique and is given 
by Eq. (^2|). For sets with dimension > 3 the bounce-free 
solution is not unique. In Eq. fl25| ) we have parametrized 
all bounce-free solutions in the 4-dimensional case. Test- 
ing the efficiency of different parameter choices on a 
model where there is only one 4 dimensional equation 
set, the s=l Heisenberg model in zero field, we find that 
even in the case where all four weights related by the 
set are equal, the most effective solution is not the most 
symmetric one (all off-diagonal a^-'s of equal magnitude). 
Thus we conclude that it is not possible to find the most 
effective solution based on the weights of a single isolated 
directed loop equation set alone. This is natural as the 
efficiency depends on the overall loop motion through all 
possible vertices and not just the motion within a certain 
subset of vertices related by the same directed loop equa- 



tion set. Our finding of the most efficient algorithm for 
the spin-1 model coincides with the most efficient direct 
algorithm found by Harada and Kawashima[^|. 

We have also shown how the directed loop equations 
can be applied to study spin-s XXZ models, lattice 
bosons, ID spinful fermions and transverse field spin 
models using the same computer code just changing the 
number of allowed states on each site as well as the ver- 
tex weights. We have also worked out expressions for 
the regions where one can construct algorithms without 
bounces for these models. 

While finding the most efficient solution to the directed 
loop equations among the many possible ones remains a 
difficult task, the solutions given here are generally more 
efficient than the solution employed in Ref . || . The many 
possible solutions to the directed loop equations should 
be seen as an asset. They are powerful tools allowing 
efficient simulations of a wide class of quantum models. 
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APPENDIX A: THE FOUR-DIMENSIONAL NO- 
BOUNCE SOLUTION WHEN TWO OR MORE 
WEIGHTS ARE EQUAL 

The no-bounce solutions in the four dimensional case 
can be parametrized as 

flu = (W 1 + W 2 -W 3 -W 4 )/2 + a 34 
ma = (Wi- W2 + W3- W 4 )/2 + a 24 (Al) 
a 23 = (-Wi + W 2 + W 3 + W 4 )/2 - (a 34 + a 24 ) 
O14 = W±— (a 34 + a 24 ) 

(A2) 

where it is assumed that W\ > W 2 > W 3 > W4. When 
two or more of these weights are equal the ordering is 
ambiguous. However, any solution obtained for a chosen 
ordering and choices of a 24 = a and a 34 = b is identical 
to a solution obtained using another ordering and other 
values of a 24 and 0,34. The purpose of Table | is to relate 
values of a 24 and a 34 for different orderings. We choose as 
a reference ordering the order ABCD which means that 
Wi = W A , W 2 = W B , W 3 = W c and W A = W D . As an 
example consider the case where Wc = Wd but Wa > 
W B > W c - Then the two orderings ABCD and ABDC 
are inequivalent. From the third entry in Table | we read 
that a choice of a 24 = a and a 34 = b for the ordering 
ABCD give the same rules as the ordering ABDC with 
the rules a 24 = Wa — a — b and a 34 = b. Wa — (—Wi + 
W 2 + W 3 + W 4 )/2. 
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TABLE I: The table shows how the solution of Eqs.flAlj) 
with ordering ABCD and 024 = a, 034 = b are related to 
other solutions with different orderings of vertices when two 
or more weights are equal. Wa = ( — Wi + W% + W3 + W4) /2. 
(* when all weights are equal the transformations (1 <-» 
2,3 <-> 4), (1 <-> 3,2 <-> 4) and (1 <-» 4,2 <-> 3) are symme- 
try transformations implying in particular that the orderings 
ABCD,BADC,CDAB and DCBA give the same rules for a 
particular choice of a and b.) 
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